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ABSTRACT 

A theory for time-dependent thermal and gas diffusion in mechanically time-rate- 
independent anisotropic poroelastic composites has been developed. This theory advances 
previous work by the latter two authors by providing for critical transverse shear through a 
three-dimensional axisymmetric formulation and using it in a new hypothesis for deter- 
mining the Biot fluid pressure- solid stress coupling factor. The derived governing equa- 
tions couple material deformation with temperature and internal pore pressure and more 
strongly couple gas diffusion and heat transfer than the previous theory. Hence the theory 
accounts for the interactions between conductive heat transfer in the porous body and con- 
vective heat carried by the mass flux through the pores. The Bubnov Galerkin finite ele- 
ment method is applied to the governing equations to transform them into a semidiscrete 
finite element system. A numerical procedure is developed to solve the coupled equations 
in the space and time domains. The method is used to simulate two high temperature tests 
involving thermal-chemical decomposition of carbon-phenolic composites. In comparison 
with measured data, the results are accurate. Moreover unlike previous work, for a single 
set of poroelastic parameters, they are consistent with two measurements in a restrained 
thermal growth test. 



Nomenclature 


Cj: degree of processing 
C p : heat capacity 
e-{. elastic strains 

e\° x \ total strains 
hg. heat enthalpy 
/zr: heat of reaction 
Kg. gas bulk modulus 

k: permeability matrix 

M : Biot’s material constant 

MWg. molecular weight 

mg. gas mass increment per unit bulk volume 

N-{. elemental shape functions 

p\ pore pressure 

q: heat flux vector 

R: universal gas constant 

T: absolute temperature 

u, w: radial and axial displacement components 
v g : average gas velocity 

oq: Biot’s material constants; pressure- stress coupling factors 

Pp thermal expansion coefficients of solid 

p g : thermal expansion coefficient of gas 

8,; unjacketed compressibilities 

(J); porosity 

k: thermal conductivity matrix 
p: gas viscosity 
p g ; gas density 

Pp ro c : density of processed solid 



p s : solid density 

Pvu-g: density of virgin solid 

X\: total stresses 

gas volume increment per unit bulk volume 
£ tot : total gas volume increment per unit bulk volume 
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1. INTRODUCTION 


Poroelasticity has been applied to numerous problems in which a fluid, diffusing 
through a deformable solid, influences the mechanical behavior in a coupled manner. The 
most common applications involving solid deformation are in geotechnical engineering 
[1]. Problems coupling fluid diffusion to thermal or electronic diffusion have also been 
solved within a rigid matrix. Such problems feature two independent field variables, 
namely temperature and pressure. Raising the level of complexity to three or more inde- 
pendent field variables, material deformation has been coupled to gas and thermal diffu- 
sion in the study of a high temperature thermal insulation material by Sullivan [2], 
Sullivan and Salamon [3], from which this work is launched, and Weiler[4], 

The theory is founded upon that of Biot [5] who developed constitutive relations 
for the elastic behavior of saturated, isotropic porous soils and Biot and Willis [6] who 
expanded them to include anisotropy. Their theory provides a mathematical model for the 
mechanical behavior of the bulk material by phenomenologically linking the interaction 
between the solid and fluid phases in the sense of "mixtures". Nur and Byerlee [7] dis- 
cussed the concept of effective stress and defined an effective stress law for isotropic, 
fluid-filled porous materials. Carroll [8] further developed the effective stress law for 
anisotropic porous materials and Kurashige [9] reported an anisotropic, thermoelastic for- 
mulation. All these formulations are for treatment of geo-materials. 

The class of problems treated here involve an evolutionary process in time during 
which solid deformation, pore pressure (due to diffusing gases) and change in temperature 
in poroelastic bodies interact. In general the process may involve fluid mass generation 
associated with conversion of a solid from a virgin state to a processed state driven by a 
change in energy level in a control volume through which mass and energy flux may trans- 
fer. In particular, the present application deals with thermally activated, chemical decom- 
position of carbon-phenolic composite material used in rocket liners. The theory is more 
general and may have other applications. 

Carbon-phenolic is a polymeric material which chemically decomposes into solid 
carbon (char) and pyrolysis gases when exposed to high temperature. Initiating at an 
exposed surface, a char layer forms, governed by chemical reaction kinetics, and advances 
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into the material lagging the thermal front. The solid carbon left behind has a porous struc- 
ture which permits generated process or decomposition gases to escape. This convective 
action assisted by a low thermal conductivity of char insulates the interior virgin material 
from heat. Hence carbon-phenolic is suitable for high temperature insulation applications. 

This paper advances the theory in [2, 3] by (1) creating a model for non uniform, 
thermal convective fields, (2) developing the coupled, three-dimensional, axisymmetric 
finite element formulation for anisotropic, axisymmetric structures in a general form and 
(3) employing a new empirical, but physical form for Biot’s pressure-stress coupling fac- 
tor. It then uses this formulation to simulate two high temperature tests done by Stokes 
[10] and, after determining a consistent set of poroelastic parameters, achieves excellent 
agreement with the test data. The following sections cover the theory and finite element 
formulation, both contain new approaches, and then apply them to experiments on decom- 
posing carbon phenolic specimens. 


2. GENERAL THEORY 

The problem is first formulated in a Cartesian frame and then set into cylindrical 
coordinates (r,9,z) such that r is the radius of a cylinder and z is the generator normal to 
the r, 0-plane where 0 is the angle of rotation about z. Axial symmetry is invoked by sub- 
jecting all fields (•) to the constraint 5(-)/00 = 0. The resultant degrees of freedom are dis- 
placements u, w in the r, z directions respectively, pore pressure p and temperature T. 

Material points are occupied by both solid and fluid, overlaid in the sense of mix- 
tures, and constitute the bulk material. The porous solid is material devoid of fluid. More- 
over solid and fluid refer to pure solid and pure fluid material, respectively. Stresses are 
usually referenced to bulk area, hence they are termed effective. However any quantity 

denoted ( ) s pertains to the pure solid; such stresses are referenced to pure solid area. Pres- 
sure p is that in the pure fluid, i.e., herein the gas. 

Mass flux through the porous solid is assumed to be slow, irrotational flow. Tem- 
perature at a material point is assumed common to both the solid and fluid, hence perfect 
heat transfer occurs between them in a control volume. Heat transfer occurs as conduction 
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in the solid and as convection by diffusion of the gas. 

The theory is formulated to achieve a linear, time-marching solution. Nonlinear 
expressions which arise are linearized. 


2.1 Constitutive Equations 

We consider a linear elastic material with porosity <|) and stiffness C whose pores 
are saturated with gas under pressure p which is taken positive for compression. The small 
elastic strain in the solid, in terms of total and thermal strains, is e x = <?j tot - Similarly 

the elastic pore fluid volume fraction change is £ = £ tot + <|)PgAr with a plus sign because 

pore fluid is lost under a temperature increase at constant pressure. The semi-complemen- 
tary form of strain energy W [11] is 

1V= ifC^-^p-ip 2 ) (1) 

from which the total stresses T, (solid plus fluid pressure referenced to the bulk area) and 
total fluid volume increment per unit bulk volume £ tot follow: 

X ‘ = l C^'-P-AD-cy, (2) 

C'“ = -^-W S AT = jjP + a^P-^T) -W S AT (3) 

where repeated subscripts are summed. Reduced indicial notation, employed to represent 
tensorial quantities, will be defined subsequently. 

In equation (3), the coefficient M is determined through the isothermal unjacketed 
test, Biot and Willis [6], and can be expressed in terms of Biot’s coefficients a, and unjack- 
eted compressibilities 8, as 

AA — ^ 8 

♦ -♦(8 1 + s 2 +8 3 )jr J +a 1 .s i jr ( (4) 
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where Kg= -V(dp/dV)j is the gas bulk modulus at constant temperature taken positive for 
decreasing volume under compression. For an ideal gas, K„~ p and the last term in equa- 
tion (3) varies as A TIT which at high temperatures will be small and may be ignored. 

The compressibilities may be determined directly from the unjacketed test, or by 
referring to Carroll [8], 




(5) 


where S s denotes an elastic compliance of the solid phase material. An alternative method 
is to first set <Xj, then use (see the dissertation by Lee [12]) 

{ 5,1 = [(",]-■( { 1 } -{<*,}) ( 6 ) 

where { 1 } is a column matrix with entries { 1 1 10 0 0 }. Conversely, given {8j}, equation 
(6) may be inverted to determine { cq } . The alternative route is adopted here because the 
freedom to determine a from experiment counters uncertainty in the values of the material 
constants during decomposition. 

For further information, details are in [12]. A general development in Cartesian 
coordinates is given by Weiler [4]; an alternative is given by Kurashige [9]. 


Axisymmetric considerations. The most general material which meets the axi- 
symmetric constraint is monoclinic with material symmetry about the r, z-plane [12] and 
is referred to here as r, z-symmetric. The elastic stiffness matrix for this material is written 
(using the notation l =rr, 2 = 00, 3 = zz, 4 = 0z, 5 = rz and 6 = r0) as 


C 11 C 12 C 13 0 C 0 

C 22 C 23 0 C 0 

C 33 0 C 0 

C44 0 c 46 

Sym C 55 0 

^66 


( 7 ) 


7 



and the material has a 4 = a 6 = 0 and p 4 = (3 6 =0. For this case, the constitutive equation (2) 
becomes 


' X 1 ' 


e\°‘ - pr ' 


a i 'i 

X 2 


4°' - pr 


tt 2 

X 3 

< 

► = [C] < 

4°' - pr 

► -p - 

«3 k 

X 4 


e tot 

e 4 


0 

X 5 


4°' - pr 


«5 

, X 6 . 


- 

0 


where 



du 

d? 



lot _ to, _ du dw 

3 dz * 5 ~ dz dr 



(9) 


2.2 Momentum Equations 


In cylindrical coordinates (r,0,z), the equations which govern the motion of the 
porous, fluid-filled material due to external loads in the absence of body forces and under 
quasi-static conditions are written as 


1 


3x i ! 9 x 6 dx 5 
dr + rdQ + dz + r {X 1 X 2 } ~ 0 


dx 


„ 1 ^ X 2 ^ X 4 1 , .. 

d7 + T W + 37 + 7 (x 6 + x 4 ) -° 


9 x 5 j 9x 4 


3x 


3 -4x e = 0 


dr r dQ dz r 5 


( 10 ) 


Axisymmetric case. Under the constraints of axisymmetric conditions (v = 0 and 
3(-)/90 = 0), £ 4 tot = e g tot = 0 from equation (9), x 4 = Xg = 0 from equation (8) and the sec- 
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ond equation in (10) is automatically satisfied. Then the equations reduce to 


a7 


^5 1 

+ 37 + 7 (t .-^=° 


dx. 3X-, 

aF + a? 


+ -x c = 0 


(ll) 


2.3 Gas Diffusion Equation 

The differential equation governing the flow of fluid through the anisotropic 
porous solid skeleton, derived from the principle of conservation of mass within a control 
volume, is 


dm g 

dt 


+ V- ( Pf F,) 


a-r 

dt 


( 12 ) 


where p g and v g are the mass density and velocity of the gas, respectively, and m g gen > 0 is 

the mass generation term (m g gen < 0 is the mass consumed). On the left hand side, the first 
term is the storage term which represents the time rate of gas mass inside the control vol- 
ume and the second term is the diffusion term which defines the flux of gas mass crossing 
the boundaries of the control volume. 

We assume the gas density p g constant during an incremental change in time. 
Hence it is treated as a parameter which is updated at each numerical time step. The incre- 
mental change in gas mass per unit bulk volume is obtained by multiplying the total vol- 
ume increment per unit bulk volume £ tot by the gas density p g . Therefore 


dm g _ ^ d? ot 
dt P 8 dt 


( 13 ) 


Hence from equation (3), the mass storage term is written in terms of pressure, 
strains and temperature as 
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dm g 

dt 


( 14 ) 


1 dp 

— ^r- +0L. 

M dt 1 


de 


tot 


dt 


dT 

-<“A + fiy-37 


In the second term on the left hand side of equation (12), the volume average gas 
velocity components in cylindrical coordinates are expressed by Darcy's law for slow, 
irrotational flow as 


v = --k • V p 

£ p 


where p is the gas viscosity and 


k = 


*1 k 6 k 5 

k 6 k 2 k 4 

k 5 k A k 3 


( 15 ) 


( 16 ) 


is the permeability matrix for an anisotropic porous body. 

The right hand side of equation (12) represents the rate of increase in gas mass per 
unit bulk volume due to some chemical conversion or phase change process and this is 
equal to the rate at which the solid phase gives up or takes in mass per unit bulk volume. 
Mathematically, 


amf" _dp 5 

dt dt 


( 17 ) 


where the solid density p s may be expressed as a linear function of the degree of process- 
ing cqas in [13], hence 


P, = 'iKirt+U-'OPproc (18 > 

in which c\ is unity at the virgin stage and zero at the completely processed stage. Expand- 
ing equation (17) using the chain rule, 

dm 8en dp s dc l dc x 

dt dc[~dT ~ ^virg Pproc)~df (19) 
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in which the rate of processing dc\/dt is defined in Section 4.1 by setting up a series of 

Arrhenius kinetic reaction equations based on experimental data. 

Substituting equation (14), (15), (16) and (19) into equation (12) gives 


Be\ ot \ 37 

-R-J-p t (“A+Wp(^) 

( 20 ) 

dc , i 

+ (p . - p ) (-T— ) - p V - (— k • V p) = 0 

v r vtrg r' proc' v df J rg r j 

where gas density is treated as a constant locally in space. Equation (20) involves temporal 
and spatial derivatives of the gas pore pressure and accounts for the changes in pore pres- 
sure due to chemical processing, temperature variations and solid deformations. 


(-=-) + p a. 


Axisymmetric case. The gas diffusion equation (20) is reduced for a monoclinic r, 
z-symmetric material in cylindrical coordinates to 


p * , d P 


M ( dt ] + P *3r 


3 r 3m 


a 


1 3r 


u 

+ a„- +a, 


dw 3w 3 m 
3 37 + “ 5 ^ 3r + 3z 

37 


dc i 

p g (ajPj + a 2 p 2 + a 3 p 3 + a 5 p 5 + <t>|3p ^ + (p virf - p proc ) w 


’ p g { rdr 


A 3/? k 5 d p 
r( TT37 + U37 ) 


44 £+££»-° 


( 21 ) 


where the coefficient M is obtained from equation (4), £4 = = 0 for r, z-symmetric mon- 


oclinic materials and equations (9) under axisymmetric conditions are employed. 


2.4 Energy Conservation Equation 


The equation of energy conservation which governs the balance of energy in a con- 
trol volume is written as 
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( 22 ) 


dE „ BE *" 1 

dt M dt 

where on the left hand side, the first term is the storage term which represents the rate at 
which energy is stored within the control volume and the second term is the diffusion term 
which defines the flux of energy convecting through the control volume and the right-hand 
side represents the rate at which energy is being generated or consumed due to the process. 

The rate of the energy stored in the control volume is the sum of the rate of the 
energy stored in the solid and gas. The first term of equation (22) may be expressed as 

dE dT 

~Si = ( 37 ) 03) 

where the subscripts s and g stand for solid and gas, respectively. 

The energy diffusion term consists of two contributions: heat conduction in the 
bulk solid and heat convection carried by the gas flow through the pores. 

*1 — flcond conv (^) 

The heat conduction term can be written by the Fourier conduction law for aniso- 
tropic materials as 

%ond = -K • VT (25) 

where K, the conductivity matrix, is similar in form to equation (16). The convective heat 
flux is the gas enthalpy /z g times the gas mass flux p g v g and is written as 

<w = V« T «-< c V f p, v , r < 26 > 

where the final term is valid if the gas is assumed to be ideal. Substituting equations (25) 
and (26) into (24) and taking the divergence of the result gives 

V-q = -V- (k- VT) -p g (C p ) VT- (Jk-Vp) - p g h g V ■ (Ik- Vp) (27) 

where the second and first forms of (26) are operated upon to get the second and third 
terms, respectively, and (15) is recalled. The reason for this juggling act is to strengthen 
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coupling between variables yet, from intuition, circumvent nonlinearities where possible. 
Notably, the last term does not appear in Sullivan and Salamon [3] because gas velocity 
was assumed to be locally constant. However for the problems involving thermal gradi- 
ents, see [ 12 ], both of these terms have a significant effect on the process. 

The energy density due to the heat generated in the control volume is expressed by 
the product of the virgin material density p^g, the heat of reaction h R> and the rate of the 

process 
dc\!dt as 


a E gen _ dc x 

~dt ~ P virg h R~^ 


(28) 


Substitution of equations (23), (27) and (28) into equation (22) yields 


[(l-MP.tcp.+tp^ty^-v. (k-vd 


1 1 rfc, 

-p (CJ VT • (-k • Vp) - ph 0 V • (-k • Vp) + p . h R - T ±= 0 

v P' q v |j, r g g rvi rg K fa 


(29) 


g ' P'g 


In the energy conservation equation (29), both pressure and temperature are inde- 
pendent variables. The nonlinear term (3rd term) can be linearized by choosing either 
pressure or temperature as a variable parameter and this choice is discussed below equa- 
tion (32). 


Axisymmetric case. The energy conservation equation (29) for the axisymmetric 
case in cylindrical coordinates is written as 


5J ( [ ( 1 - ♦> P, (C„), + *p, (C„) J T,) + p 


dc. 


P' s 


p' g J ‘ g’ 


virg 


l R 


dt 


1_9 

r dr 


dT dT 1 

r{K 'Tr +K ^ ) 


3 , 3T 3T 

Tz (K S3? + K iTz ) 


-P,< C ,) S 

, , 1 d 
- p A <737 


k\dp k$ dp dT £5 dp k^dp dT 
|l dr + |l dz dr + jj. dr + p. dz ^ dz J 


k\ dp k 5 dp 
r ( — -=r- + — ^-) 
p. dr p .dz 


d k 5 dp k-i dp 

+ t-(-/ + -^-)} =0 

dz p. dr p.dz 


( 30 ) 


13 



Because this equation employs the final term in equation (26), it is restricted by the 
assumption of an ideal gas. 


3. AXISYMMETRIC FINITE ELEMENT FORMULATION 


The axisymmetric finite element equations are obtained by successively applying 
the Galerkin method to equations (11), (21) and (30) in the following form: 


J/V- (LHS_of_referenced_equation) rdA = 0 (31) 

A 

where LHS denotes the left hand side and A is the domain in the r, z plane. 

The nonlinear, penultimate term in equation (30), namely, 


-p ^ c p \ 


k\ dp £5 3d dT &5 dp k$dp dT~ 

|X dr + JJ. dz dr + (J. dr + |l dz dz _ 


(32) 


is linearized by choosing the spatial derivatives of temperature as parameters since pres- 
sure is a more sensitive factor to the processes studied. Then applying Green's theorem, 
and integrating (31) and (32) by parts (Lee, 1993), the following matrix equation is 
obtained. 


[C]j t {a] + [K] { a } = { F } 


where 


[ C ] 


0 0 0 0 

0 0 0 0 

C pu C pw C pp C pT 

0 0 0 Cjj 

W = (>. K Fp F r] 


[*•] = 


^uu K uw K up K u j 
^ wu L K W p ^ wT\ 


0 


0 Knn 0 

PP 


0 0 K Jp Kjj 


{ a } - [u w p i] 


The elements of [C], [AT] and {F} are explicitly written in the Appendix. 


(33) 


(34) 
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It should be noted that the coupling term Kj p does not appear in Sullivan and Sala- 
mon [3]. They treat the spatial derivatives of pressure as parameters which causes these 
terms to be deposited on the diagonal of the [K] matrix rather than in the K Tp term. The 

above formulation results in a more strongly coupled system of equations and is employed 
here. 

For integration over time, the variably weighted Euler numerical method [14] is 
employed which leads to 




(35) 


where 


_ 1 _ 

At ' 1 


[K eff V = —JC] n + Q[K] n -Q[H] n 


(36) 


{ F eff } n = { F} n + 


1 

At" 


[C] n - (1-9) [K] n -Q[H] n 


{a}" 


(37) 


and superscript n denotes the time step. In equations (36) and (37), the matrix [H] may be 
defined as 


3{F) n 

( 38 ) 

d{a} 

For solution by the fully implicit time integration scheme, the matrix [H] is required. 
Details of its formulation are given in the dissertation by Lee [12]. Nonzero elements of 
[H] are given in the Appendix. 

In equation (35), [Af e ff] n and [F e g-} n are the effective stiffness matrix and effective 
force vector at time step t n respectively. To obtain the material displacements, pore pres- 
sure and temperature at each time step, equation (35) is solved for [a} n+1 simultaneously 
using a factorization procedure which is a version of Gauss elimination with partial pivot- 
ing [15], 


15 



4. DECOMPOSING CARBON PHENOLICS 


The governing equations for poroelastic material with thermal and gas diffusion 
provided in Section 2 are applicable to thermochemically decomposing carbon-phenolic. 
The momentum equations remain the same, but the generation terms in the gas diffusion 
and the energy conservation equations require further specialization. The objective is to 
specialize and apply the linearized theory to models of carbon-phenolic material and sim- 
ulate laboratory tests in order to demonstrate it and establish poroelastic parameters over a 
nonlinear decomposition process. 


4.1 Governing Equations 


For the mass generation term (19) in the gas diffusion equation (12), the rate of gas 
mass accumulation is due to decomposition reactions and, from Sullivan and Salamon [3], 
is given by 


a-r 

dt 


N 


i = 1 


dic l ) i 

dt 


( 39 ) 


where RF is the weight fraction of the virgin composite which is resin, p virg is the virgin 
density of the composite, (W r 0 ) i is the fractional weight of the resin which undergoes the 

1 th reaction, (W c )i is the fractional weight of the resin which is left as a solid residue by the 

i -th reaction, and N is the number of reactions. Following [3], the rate of charring d(c\\/dt 
is expressed by the Arrhenius kinetic reaction equation 


die i). 
dt 





( 40 ) 


where A s and n s are the Arrhenius constants and E s a is the activation energy for the chem- 
ical reaction. These constants for carbon-phenolic are listed in Table 1 . 
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In this work, pyrolysis gases are treated as ideal and the gas density p g in the gov- 
erning equations is obtained by the ideal gas law 

MW 

Pg rT p (41) 

where MWg is the molecular weight of gas and R is the universal gas constant. The 
mechanical and thermal properties of dry carbon-phenolic are listed in Table 2. 


4.2 RTG and FTE Tests and Their Finite Element Models 

Stokes [10] conducted two high temperature experiments on cylindrical carbon- 
phenolic specimens (dry FM5055), fabricated so that the plane of the carbon fabric is per- 
pendicular to the longitudinal axis of the specimen, in order to determine thermomechani- 
cal behavior during chemical decomposition. The restrained thermal growth (RTG) test 
measured both the stress required to hold the specimens at a constant longitudinal strain 
and the resulting lateral strain. The free thermal expansion (FTE) test measured the result- 
ing longitudinal strain. The specimens were heated uniformly at a rate of 5.55 degrees 
Kelvin per second (10 degrees Fahrenheit per second), and this was controlled using ther- 
mocouples embedded in the specimens. 

The geometry, coordinate system and finite element mesh for the RTG and FTE 
test specimens are shown in Figure 1. In the RTG test, it was assumed that the end con- 
straints prevent pyrolysis gases from flowing along the axial z-direction. Hence, imperme- 
able conditions dp/dz = 0 are prescribed at z = 0 and z = b. However in the FTE test, the 
pyrolysis gases do flow along the z-direction, hence to allow the gases to escape, atmo- 
spheric pressure is specified at z = 0 and z = b. The initial conditions are T = 293 degrees 
Kelvin and p = 1 atm in the model and the boundary conditions are summarized in Table 3. 

The energy generation term (28) in equation (22) is simplified in order to maintain 
a uniform temperature increase in the model. An internal heat source is prescribed in each 
finite element through the integral equation 
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dE gen 

dt 


(QinJi = \ N i<lint rdA 

A 


( 42 ) 


where £?j nt is the heat generation per unit area, and (£?int)i ‘ s an equivalent nodal value. The 
value of 9 in t , found by trial and error in order to achieve the rate of temperature increase of 

the specimen (5.55 degrees Kelvin per second), was determined to be 11 .4 x Joule! 
rrt^sec. 


4.3 Results for the pressure-stress coupling factors a 

Calculation of the pressure-stress coupling factors from equations (5) and (6) for 
the virgin elastic composite material requires an experiment, the Biot and Willis [5, 6] 
unjacketed test. Unfortunately such an experiment for carbon-phenolic has not been con- 
ducted. Moreover after decomposition commences, (5) and (6) are suspect because of 
severe changes in material morphology and properties. Hence further experiments would 
be desireable. In short, requisite data to calculate these poroelastic parameters does not 
exist 

Results are obtained by simulating the RTG and FTE tests and varying a, then 
comparing the material response to experimental data. This was done by physical intuition 
and trial and error based upon a hypothesis introduced by Lee [12]. This hypothesis 
assumes (1) isotropic porosity, i.e., oq = a, i = 1, 2, 3, 5, (2) a higher proportion of closed 

pores in virgin material, and (3) increased pore channel opening under high shearing stress 
action. Consequently, the pressure-stress coupling factor (1) is initially high and, if data is 
available, can be determined from material properties and equations (5) and (6), (2) 
decreases as porosity opens, thus inversely as permeability increases, and (3) is further 
decreased if shear stress enhances pore channel growth. (Note: the procedure employed 
here is to guess a, then use equations (6) and (4) to compute M; hence a is a curve fit 
parameter.) 

In Figures 2 and 3, average restraining stress and average lateral strain of the RTG 
test simulation are plotted versus average temperature. Also, in Figure 4, average longitu- 
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dinal strain of the FTE test simulation is plotted versus average temperature. The variation 
with temperature of the pressure-stress coupling factor a used in the RTG and FTE simu- 
lations, shown in Figure 5, is that for which the numerical solutions most closely match 
the measured responses. Only values at comers of the curves were adjusted, those in 
between corner values were linearly interpolated. The pressure-stress coupling factor for 
the RTG test simulation was decreased more than that for the FTE test simulation due to a 
higher shear stress level in the RTG test simulation. Notably for the RTG test simulation, a 
single choice for the varying pressure-stress coupling factor correlates with data for both 
restraining stress and lateral strain; Sullivan and Salamon [3] required different constant 
values to achieve satisfactory results in their parametric study. 


5. CONCLUSION 

A coupled set of governing equations for a poroelastic solid with thermal and mass 
diffusion is derived in a practicable form which requires a minimum of experimental data, 
yet are sufficiently realistic to entertain engineering problems. In the energy conservation 
equation, both conductive heat transfer in the bulk solid and convective heat transfer car- 
ried by the diffusive gas through the bulk porous solid are treated. The convection terms in 
which pressure is the independent variable are included in the energy conservation equa- 
tion so that the convection heat transfer by the diffusive gas is accounted for more accu- 
rately than in [3]. These convection terms are coupled with the gas mass diffusion term in 
the gas diffusion equation, and therefore pressure appears in both the gas diffusion equa- 
tion and the energy conservation equation as an independent variable. Hence the theory 
now handles the thermal and gas diffusion in the poroelastic solid in a strongly coupled 
manner. 

The material formulation includes the most general anisotropic material under the 
constraints of axisymmetry, and therefore permits application to realistic, three-dimen- 
sional composite structures. Importantly the axisymmetric formulation permits accurate 
computation of the transverse shear stress which is hypothesized to play a significant role 
in pore channel opening during material decomposition and in turn enables finer definition 
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of the pressure-stress coupling factor during decomposition. 

The coupled theory, when applied to carbon-phenolic material tests, provides con- 
sistent and close correlations with the experimental data using an empirically determined 
pressure-stress coupling factor suggested by Lee [12]. The improved correlation with 
experimental data in simulating these tests is attributed to the new formulation for the 
pressure-stress coupling factor. 

It is anticipated that the strong coupling of the equations will play a major role in 
accurate computation of spatially nonuniform thermal problems such as those encountered 
in rocket liner structures during firing. 
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APPENDIX 

The elements of [C], [A] and {F} are 


(*■„„) ij = \N,.r (C u Nj, r + C u iN, + C 15 NjJ rdA 

A 

+ i" i .ACu"l,, + CvTNj+C sa Nj',) rdA 

A 

+ f»i ^n N i,r + C li \N j + C^N jt z ) dA 
A 

(*„,,) , = 1 W,. , (C I5 W ; , , + C 13 W. : ,) rdA + J/V t , (C 55 W ;> , + C 35 W ; . ,) rdA 

A A 

+ jw j (C 25 N y „ + C 23 W /I )dA 

A 

(*„,)« = "/dA -jA^a^rdA -Jw.a/TdA 

AAA 

<*.!•)(■ = -K r e 1 A'/dA-Jw j , z e 5 w / dA-|/v.e 2 w 3 dA 

AAA 

(*»„) ,J = K r (C,sty, , + ClS^j + CjjAO. ,) "M 

A 

+ K, (C i 3 Hj, r + C a -Nj + c 35 /v. ,) rdA 

A 

= Kr<C 55 W y ,, + C 33 ^ ! )rdA+ R^AT, + C 33 Ar ; )rdA 
A A 

(*„„>„ = -Jw i , r a 5 W / dA-Jw,. 2 a 3 W / dA 

A A 

(K wT hi = -K,e 5 ^A-|)V, ! G3W / dA 

A A 

k k k k 

(K ) .. = f/V. n (— N ir +—N iz )rdA+ ftf, zP<j (— M r + — AF 7 ) rdA 
v PP' i; J l > rr g (l h r [i J' z J i » 2K <? v (l J' r fi /> z ' 


(A-l) 


(A-2) 

(A-3) 

(A-4) 

(A-5) 

(A-6) 

(A-7) 

(A-8) 

(A-9) 
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(A- 10) 



(K u ) = Jn. r ( K^. r + k s N jt 2 ) rdA + jN ltZ (k 5 /V._ r + k z ) rdA 

A A 

(C P« } w = JXP* ( a i N j, r + a 5 N J, z) rdA + \ N P 2 P g N j dA 
A A 

(C„J = l N fi^ a 5 N i,r + ^ jf! )rdA 


< C PP ),7 = fa&fU 

A 

(Cpj) ■ • = ( a iPl ^2^2 ^3^3 *** ^5^5 N^rdA 

A 

< c rr), y = \NiW--*)V,(C t ) +*>,<£,) JNfdA 

A 

(f«), = 

(f»), - Jw,v* 

s » 

( p p ) i = - (Pv,>,“PprjJ N i^l) rdA -pJ^i( V I)* n ik r ^ 

A i p 

( p rh = " Pv«> s M N i e i rdA ~ 

A s y 

where 


(A-ll) 
(A- 12) 
(A- 13) 

(A- 14) 
(A- 15) 
(A- 16) 
(A- 17) 
(A-18) 
(A-19) 
(A-20) 


Q\ - ^llPi + ^12^2 + ^13^3 + ^15^5 
2 3 = C B p 1 + C 23 ( 3 2 + C 33 P 3 + C 35 P 5 


Q 2 — ^12^1 + ^22^2 + ^23^3 + ^25^5 
Qs = C 15 Pj + C 2 5p 2 + ^35(^3 + ^55(^5 


and where v g * 


is the volume average gas velocity specified along the external boundary of 
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arc length s p , q* is the heat flux specified on s T , n x represents components of the outward 

unit vector normal to the boundaries, and t x and t z represent traction components distrib- 
uted over the external boundary of arc length s u and s w in the r and z directions, respec- 
tively. 

The [H] matrix is written explicitly as 


where 


[H] 


0 0 0 0 
0 0 0 0 
0 0 0 H pT 
000 % 


(A-22) 


A 

(A-23) 

=-p 

(A-24) 
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Table 1: Constants in the decomposition model for dry carbon-phenolic. 


Reaction 
number (i) 

E s a (J /mole) 

A s (1/sec) 

n s 

^0 


1 

88764.4 

1.207305 x 10 10 

3.5 

0.0015 

0 

2 

117236 

4.057500 x 10 9 

6.5 

0.095 

0 

3 

211443.5 

3.857777 x 10 14 

6.5 

0.59 

0.29 

4 

272155 

5.583611 x 10 15 

3.3 

0.3 

0.19 


RF = 0.3, h R = 0, MW g = 0.03 kg/mole 
































Table 2 : Mechanical and thermal properties of dry carbon-phenolic. 


property 

temperature range 

T < 450 °K 

450 °K < T < 533 °K 

T > 533 °K 

e l 

1.517X10 10 Pa 

-1.412x10 s T + 7.87X10 10 Pa 

3.447X10 9 Pa 

e t 

1.793X10 10 Pa 

-1.7449x10 s T +9.645xl0 10 
Pa 

3.447xl0 9 Pa 

V TL 

0 

-0.00343 T+ 1.83819 

0 

V LT 

0 

-0.00289 T + 1.551 

0 

Vj 

0 

0 

0 


property 

value 

property 

value 

Pt 

0.000006 m/m-°K 

bL 

0.000012 m/m-°K 

(^T)virg /M* 

5xl0‘ 21 m 3 -sec/kg 

(^L)virg /M- 

5x1 O' 21 m 3 - sec/kg 

(Mchar /M- 

5xl0' 13 m 3 -sec/kg 

(^L)char /P- 

5x1 0' 1 3 m 3 - sec/kg 

^virg 

0.02 

^char 

0.20 

Pvirg 

1500.0 kg/m 

Pchar 

1300.0 kg/m 3 

(Cp) solid 

1400.0 J/kg-°K 

(Cp)gas 

1088.0 J/kg-°K 

K soiid 

1.44 J/m-sec-°K 

K gas 

0.0 J/m-sec-°K 
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Table 3 : Boundary conditions for the RTG and FTE test simulations. 


location 

RTG Test 

FTE Test 

r = 0 

3r dp 

U = = 0 

dr dr 

dT dp 

“ dr ~ dr ° 

n 

dT „ 

= 0, p = 1 atm 
dr 

3 T 

-*r- = 0, p = 1 atm 
dr 

1 = 0 

dT dp 

w = 3z = 3. = 0 

dT n 

w = = 0, p = 1 atm 

dz 

— 

H 

dT dp 

W = 3z 3z ° 

dT r, 

^ = 0,p = 1 atm 
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Dimension (cm) 



a 

a 

RTG 

0.6350 

2.54 



2.54 






Figure 1 Finite element mesh for the RTG and FTE test simulations. 
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Figure 2 Restraining stress versus temperature in the RTG test simulation for 
varying pressure-stress coupling factor. 
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Figure 3 Lateral strain versus temperature in the RTG test simulation for 
varying pressure-stress coupling factor. 
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Figure 4 Longitudinal strain versus temperature in the FTE test simulation 
for varying pressure-stress coupling factor. 



Pressure-Stress Coupling Factor 



Figure 5 Typical variation of the pressure-stress coupling factor with 
temperature for the RTG and FTE test simulations. 
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